close all; clc; clear;
% This file generates the impulse response presented in Figure 1 and
% Figure 2 in Section 2 of the paper. 

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Parameters and Steady-State Values %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


beta                    =   0.99;   % Discount factor.
chi                     =   0;      % Inverse of Frisch elasticity. 
sigma                   =   1;      % Inverse of EIS.

phi                     =   0.75;   % Prob. firms cannot adjust prices. 

rho                     =   0.9;    % Persistence of natural rate shock. 

% Slope Phillips Curve:
gamma                   =   ((1-phi)*(1-phi*beta)/phi)*(sigma+chi);

T                       =   21+1;   % Number of periods.           

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%             Creating IRFs                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

H                       =   3;      % Peg length. 

% Pre-allocating:
pi0                     =   zeros(T,1);
pi1                     =   zeros(T,1);
x0                      =   zeros(T,1);
x1                      =   zeros(T,1);
int                     =   zeros(T,1);
int1                    =   zeros(T,1);
a                       =   zeros(T,1);

% Initial values:
int(1:H)                =   -(1/beta-1)*ones(H,1);
int1(1:H)               =   -(1/beta-1)*ones(H,1);
x0(H)                   =   -(1/sigma)*int(H);
pi0(H)                  =   gamma*x0(H);
pi1(H)                  =   gamma*x1(H);

for jx = 1:H-1
    x0(H-jx)                =   x0(H-jx+1) - (1/sigma)*int(H-jx) + (1/sigma)*pi0(H-jx+1);
    pi0(H-jx)               = gamma*x0(H-jx) + beta*pi0(H-jx+1);
end
y0                      =   x0;

% Now hit it with a one unit shock to a(t)
for jx = 1:T
    a(jx,1)                 =   rho^(jx-1)*1;
end
yf                      =   ((1+chi)/(sigma + chi))*a;
rf                      =   (sigma*(1+chi)*(rho-1)/(sigma + chi))*a;
x1(H)                   =   -(1/sigma)*int1(H) + (1/sigma)*rf(H);
int1(H+1:end)           =   rf(H+1:end);

% Creating IRFs:
for jx = 1:H-1
    x1(H-jx)                =   x1(H-jx+1) - (1/sigma)*int1(H-jx)...
                                + (1/sigma)*pi1(H-jx+1) + (1/sigma)*rf(H-jx);
    pi1(H-jx)               =   gamma*x1(H-jx) + beta*pi1(H-jx+1);
end

y1                      =   x1 + yf;
irfy3                   =   y1-y0;
irfpi3                  =   pi1-pi0;

%% Different Peg. 
H                       =   6;

% Pre-allocating:
pi0                     =   zeros(T,1);
pi1                     =   zeros(T,1);
x0                      =   zeros(T,1);
x1                      =   zeros(T,1);
int                     =   zeros(T,1);
int1                    =   zeros(T,1);
a                       =   zeros(T,1);

% Initial values:
int(1:H)                =   -(1/beta-1)*ones(H,1);
int1(1:H)               =   -(1/beta-1)*ones(H,1);
x0(H)                   =   -(1/sigma)*int(H);
pi0(H)                  =   gamma*x0(H);
pi1(H)                  =   gamma*x1(H);

for jx = 1:H-1
    x0(H-jx)                =   x0(H-jx+1) - (1/sigma)*int(H-jx) + (1/sigma)*pi0(H-jx+1);
    pi0(H-jx)               = gamma*x0(H-jx) + beta*pi0(H-jx+1);
end
y0                      =   x0;

% Now hit it with a one unit shock to a(t)
for jx = 1:T
    a(jx,1)                 =   rho^(jx-1)*1;
end
yf                      =   ((1+chi)/(sigma + chi))*a;
rf                      =   (sigma*(1+chi)*(rho-1)/(sigma + chi))*a;
x1(H)                   =   -(1/sigma)*int1(H) + (1/sigma)*rf(H);
int1(H+1:end)           =   rf(H+1:end);

% Creating IRFs:
for jx = 1:H-1
    x1(H-jx)                =   x1(H-jx+1) - (1/sigma)*int1(H-jx)...
                                + (1/sigma)*pi1(H-jx+1) + (1/sigma)*rf(H-jx);
    pi1(H-jx)               =   gamma*x1(H-jx) + beta*pi1(H-jx+1);
end


y1                      =   x1 + yf;
irfy6                   =   y1-y0;
irfpi6                  =   pi1-pi0;

%% Different Peg. 
H                       =   10;

% Pre-allocating:
pi0                     =   zeros(T,1);
pi1                     =   zeros(T,1);
x0                      =   zeros(T,1);
x1                      =   zeros(T,1);
int                     =   zeros(T,1);
int1                    =   zeros(T,1);
a                       =   zeros(T,1);

% Initial values:
int(1:H)                =   -(1/beta-1)*ones(H,1);
int1(1:H)               =   -(1/beta-1)*ones(H,1);
x0(H)                   =   -(1/sigma)*int(H);
pi0(H)                  =   gamma*x0(H);
pi1(H)                  =   gamma*x1(H);

for jx = 1:H-1
    x0(H-jx)                =   x0(H-jx+1) - (1/sigma)*int(H-jx) + (1/sigma)*pi0(H-jx+1);
    pi0(H-jx)               = gamma*x0(H-jx) + beta*pi0(H-jx+1);
end
y0                      =   x0;

% Now hit it with a one unit shock to a(t)
for jx = 1:T
    a(jx,1)                 =   rho^(jx-1)*1;
end
yf                      =   ((1+chi)/(sigma + chi))*a;
rf                      =   (sigma*(1+chi)*(rho-1)/(sigma + chi))*a;
x1(H)                   =   -(1/sigma)*int1(H) + (1/sigma)*rf(H);
int1(H+1:end)           =   rf(H+1:end);

% Creating IRFs:
for jx = 1:H-1
    x1(H-jx)                =   x1(H-jx+1) - (1/sigma)*int1(H-jx)...
                                + (1/sigma)*pi1(H-jx+1) + (1/sigma)*rf(H-jx);
    pi1(H-jx)               =   gamma*x1(H-jx) + beta*pi1(H-jx+1);
end


y1                      =   x1 + yf;
irfy10                  =   y1-y0;
irfpi10                 =   pi1-pi0;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%           Plotting results                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t                       =   0:T-2;

pif                     =   zeros(T,1);

figure
plot(t,yf(1:T-1),'-b',t,irfy3(1:T-1),'--b',t,irfy6(1:T-1),':b',t,irfy10(1:T-1),'-.b','Linewidth',1.5)
xlabel('Horizon','Interpreter','latex','fontSize',14)
ylabel('Output','Interpreter','latex','fontSize',14)
h = legend('H=0','H=3','H=6','H=10','Location','SouthEast');
set(h,'Interpreter','latex')
grid on
set(gca,'layer','top')

figure 
plot(t,pif(2:T),'-b',t,irfpi3(2:T),'--b',t,irfpi6(2:T),':b',t,irfpi10(2:T),'-.b','Linewidth',1.5)
xlabel('Horizon','Interpreter','latex','fontSize',14)
axis([0 T-2 -0.5 0.1 ])
ylabel('Expected Inflation','Interpreter','latex','fontSize',14)
h = legend('H=0','H=3','H=6','H=10','Location','SouthEast');
set(h,'Interpreter','latex')
grid on
set(gca,'layer','top')


